using Pkg
Pkg.develop(path = "~/Documents/EcologicalNetworksDynamics.jl-main/")This section of the tutorials initiates you into using the Bioenergetic Food Web Model.
Getting the Julia Package
Currently Not Available for Guest Users
Navigate to the GitHub location for the BioEnergetic Model.
It should look like this:
Then download to your Documents Folder the package.
Once this folder is unzipped in your Documents Folder, you can install it using Pkg.develop.
And now, you are ready to construct a script using the modelling tools!
## My first BEFW Modelling
## Packages I need
using DataFrames, Plots, Random, Distributions
using EcologicalNetworks, EcologicalNetworksPlots, EcologicalNetworksDynamics
## Time to do some Experiments!Using the Julia Package
Preamble: The Bioenergetic Food Web.
It is very worth pausing for a moment and visting this paper by Dr. Eva Delmas. It describes the BioEnergetic Food Web Model, provides some history of the model, and also showcases how the original Julia version of the model worked. This section of the tutorials is designed to introduce a newer, faster and more flexible version of the Julia package.
A very basic interpretation of the model is as follows:
- The model is a model of biomass dynamics, not numbers of indiduals.
- The model is comprised of an equation for plants (producers) an equation for consumers (herbivores, predators).
- Plants have traits that make their biomass grow and sets their carrying capacity; they are eaten by consumers via an equation describing a functional response. To link to ecology courses you have, this can be logistic growth for the plant and a type II functional response.
- Consumer have three sets of traits. One is metabolism, which is a rate that describes losses of biomass due to, well, metabolism! The second set of traits correspond to the functional response - for example describing attack rates and handling times of prey items. The third corresponds to the numerical response, or the conversion of biomass they eat into new biomass (e.g. babies)
- We can make complex networks and systems of these equations by letting many of these paramters scale with body size via the rules of allometry and the Metabolic Theory of Ecology. This trick expands the two equations to n = number of species when we provide the toolbox with a list of species biomasses.
- Embedded in this process are rules about how to distribute species of different sizes across trophic levels, so that we end up with predators, omnivores, herbivores and plants. We can also designate the body size relationships for different groups of organisms, like producers, invertebrates, endothermic vertebrates and ectothermic vertebrates.
- Once we’ve done this, we can simulate the biomass dynamics of complext communities. And we can summarise things like changes in biodiversity (number of species), stability (coefficient of variation of the time series) and anything about the biomass of species, trophic levels or the community!
Preamble: Setup
One of main advantages of running food web models in Julia is that simulations are fast and can be readily stored in your active project. With this in mind, make a new folder in your project called out_objects (right click > New Folder).
A first run of the Ecological Network Dynamics model (BEFW)
There are four major steps when running the BioEnergetic Food Web model in Julia. These should be familiar from our introduction to the DifferentialEquations package:
- Generate an initial food web network
- Set the parameters for each species in the network to generate the equations
- Simulate the network and equations
- Explore output and plot
While in the previous example with Differential Equations we assumed a simple 2-species network, one of the new activities here is to take advantage of a rich history of theory and tools to construct species rich networks with appropriate structural properties, such as connectance/complexity and levels of generalism/specialism and things the number of trophic levels and a body size distribtion of the species across trophic levels.
Step 1: generate an initial network
Here we make a foodweb, actually, a food chain, from an adjacency matrix with the FoodWeb method.
A = [0 0 0; 1 0 0; 0 1 0] # 1 basal producer ⋅ Species 2 eats 1 ⋅ Species 3 eats 2
foodweb = FoodWeb(A)FoodWeb of 3 species:
A: sparse matrix with 2 links
M: [1.0, 1.0, 1.0]
metabolic_class: 1 producers, 2 invertebrates, 0 vertebrates
method: unspecified
species: [s1, s2, s3]Step 2: Generate the model parameters
Once the foodweb is created, the next step is to attribute values to the model parameters. This can be simply done by calling the method ModelParameters with foodweb as an argument.
# construct the equations and fixed parameters
# see below for body size dependent parameters etc
params = ModelParameters(foodweb)ModelParameters{BioenergeticResponse}:
network: FoodWeb(S=3, L=2)
environment: Environment(K=[1, nothing, nothing], T=293.15K)
biorates: BioRates(d, r, x, y, e)
functional_response: BioenergeticResponse
producer_competition: ProducerCompetition((3, 3) matrix)
temperature_response: NoTemperatureResponseStep 3: Simulate biomass dynamics
Everything is ready to run the simulation, which can be simply done by calling simulate with the model parameters (params) and a vector species’ initial biomass (B0).
# create body sizes for each species
B0 = [0.5, 0.5, 0.5]
# use simulate function
# builds equations and uses DiffEq to run them!
sim = simulate(params, B0)retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 129-element Vector{Float64}:
0.0
0.07219249113352678
0.18654819333861794
0.31356331425611517
0.47345866183962954
0.6641617388942105
0.9131302533115802
1.2562392740181836
1.6309481665915588
2.089143773957029
2.636128104609133
3.293883734334714
4.047276405948346
⋮
452.3138015430547
456.8509459661425
461.5266562665336
466.1513481974734
470.70605617505174
475.3355823765842
480.0140185879458
484.5101300067372
489.02575718421326
493.69451113794287
498.331767667838
500.0
u: 129-element Vector{Vector{Float64}}:
[0.5, 0.5, 0.5]
[0.4272316681574567, 0.47585426785094076, 0.534072293177502]
[0.3473341401191793, 0.42303815072107315, 0.5859773803469567]
[0.294791058662908, 0.35687684744418796, 0.635235407219315]
[0.26038496182613113, 0.28139632314454743, 0.6781234242711596]
[0.24356863933697817, 0.2144660212785269, 0.7011935826677415]
[0.24075042789862067, 0.15961205273400708, 0.699135084180834]
[0.25337593682400383, 0.11902051844449948, 0.6668360012114383]
[0.2774846493473314, 0.097203714412661, 0.6177409318953245]
[0.31307920482758345, 0.0857704806321114, 0.5548733982179885]
[0.3554414132717484, 0.0839628299185932, 0.48528724408580537]
[0.39365248269575775, 0.09234997261222119, 0.4145536079678415]
[0.40739756802826455, 0.11063681962703542, 0.35246966612900116]
⋮
[0.22059513756843632, 0.18893233853869676, 0.04886706825052883]
[0.22055285427385707, 0.1890161486143565, 0.04885919401739655]
[0.22057127277943078, 0.18907642814450482, 0.04884560373332121]
[0.22065930151320282, 0.1890720055670567, 0.048832142934495525]
[0.22074484049723223, 0.18900107865719, 0.04883177554012091]
[0.22078067645052313, 0.18891713876325852, 0.04884189574210531]
[0.22073129917064033, 0.18885864940545732, 0.048859704977389826]
[0.2206062055609253, 0.18892017802287123, 0.048867499198903715]
[0.22055510500794875, 0.18900674970792372, 0.04886045346488585]
[0.2205647417518207, 0.18907108625114374, 0.048847621938152584]
[0.22064629570400218, 0.18907728258234632, 0.04883326864755859]
[0.22059687801237501, 0.1890242584230083, 0.04885048633720308]Step 4: Seeing the outputs!
Eventually you may want to plot the biomass dynamics - the trajectory - of your community to see what is happening. For our minimal example, it can be done as follows:
# create multiple objects: time = t pulled frpom the sim.t component
# and Bx = biomass for each species pulled from the larger sim object
# note how julia allows multiple things on left of the =
t, B1, B2, B3 = sim.t, sim[1,:], sim[2,:], sim[3,:]; # unpack variables
# Plot the basal species
plot(t, B1, lw = 3, label="Producer", xlabel = "Time", ylabel = "Biomass")
# add the herbivore
plot!(t, B2, lw = 3, label="Consumer")
# add the top predator
plot!(t, B3, lw = 3, label="Top consumer")
More generally, you could also just plot the actual simulate object directly. Much easier!
plot(sim, label = ["Producer" "Consumer" "Top consumer";])
A More Complex Example
Step 1: Generate the initial network
In order to run the BEFW model with a more complex network, we have to construct an initial food web network (an adjacency matrix) using the niche model. The network is characterised by the number of species in the network and its connectance/complexity value.
Note that we are now using functionality provided by the EcologicalNetworks package.
S = 20; # define the number of species
C = 0.2; # define the connectance (complexity) of the network
# construct the food web
Random.seed!(12325) # ensures your network and this one are the same
foodweb_niche = FoodWeb(nichemodel, S, C=C)
# see it:
foodweb_niche.A20×20 SparseArrays.SparseMatrixCSC{Bool, Int64} with 80 stored entries:
⢠⡀⠀⠀⠀⠀⠀⠀⠀⠀
⢐⣊⣀⠀⠀⠀⠀⠀⠀⠀
⠴⣶⣶⣤⣀⠀⠀⠀⠀⠀
⠈⠩⠿⣿⠄⠀⠀⠀⠀⠀
⠉⠩⠭⠭⣟⡓⠒⠒⠂⠀Step 2. Setting up the paramters, body masses (species) and running the model!
As above, our next step is to define a vector of bodymasses and then pass this, and the network to the simulate function. Here we combine the Uniform function from the Distributions package with the rand function from the Random package.
# construct the equations and fixed parameters
# see below for body size dependent parameters etc
params_niche = ModelParameters(foodweb_niche)
# define bodymasses between 0 and 1 and get S = 20 of them.
# this will ensure your plot looks like the one in the document
Random.seed!(123)
B0 = rand(S)
# simulate using params and bodymasses
# note additional argument tmax for max time steps
# default is 300
sim_niche = simulate(params_niche, B0)┌ Info: Species [5] went extinct at time t = 33.834970465767604.
└ 1 out of 20 species are extinct.
┌ Info: Species [10, 8] went extinct at time t = 49.801980626257766.
└ 3 out of 20 species are extinct.
┌ Info: Species [13] went extinct at time t = 67.54842232722383.
└ 4 out of 20 species are extinct.
┌ Info: Species [14] went extinct at time t = 76.05686387794776.
└ 5 out of 20 species are extinct.
┌ Info: Species [12] went extinct at time t = 94.25974163851537.
└ 6 out of 20 species are extinct.
┌ Info: Species [15, 19] went extinct at time t = 180.6781729701157.
└ 8 out of 20 species are extinct.retcode: Terminated
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 58-element Vector{Float64}:
0.0
0.05642845716276844
0.11891649263385903
0.20769683238716719
0.3045875314442914
0.44216128503645374
0.6112783739340873
0.8303542080772176
1.0917849224507885
1.4340233039417565
1.8447846396665963
2.3653444519946065
3.0115146437752442
⋮
103.19272243700584
111.88147068621589
122.54731908523809
133.74922122845595
147.26632373648013
162.38512727362877
180.6781729701157
180.6781729701157
201.57862868705948
226.71369483198157
254.89823928856072
289.09217901466906
u: 58-element Vector{Vector{Float64}}:
[0.906299638797481, 0.44349373245960455, 0.7456733811393941, 0.5120830400366143, 0.2538490889415096, 0.33415153638191886, 0.4273278808735992, 0.867547200255958, 0.09913361484360417, 0.12528740769155033, 0.6922086620547391, 0.13655147513745736, 0.03209667335274724, 0.3505458214588266, 0.9303323763821093, 0.9594335994071538, 0.5819123423876457, 0.3114475007050529, 0.12114752051812694, 0.20452981732035946]
[0.8506750248039066, 0.316016267593785, 0.651840748854925, 0.4766336454051591, 0.25339023877433103, 0.3308187673571854, 0.3608638414554032, 0.7307871678092578, 0.10181476103833324, 0.1302817851708985, 0.7090830957307426, 0.143047476716464, 0.033895623646412715, 0.3678562194543374, 0.9644776003724441, 1.0175145940496193, 0.6126719605021673, 0.3314093482768925, 0.12683785756057048, 0.2129203618593246]
[0.79061537214418, 0.23954384290971395, 0.5567140397232614, 0.4344170769662325, 0.24986737440786716, 0.32063982029985305, 0.2984775562116461, 0.5956422271584304, 0.10443600830700829, 0.1346128355925403, 0.7185118003079496, 0.14901914689388757, 0.03558419101698233, 0.38339997438386947, 0.9929266560084749, 1.06953278911396, 0.6433589739159863, 0.3561234135700153, 0.1320907952178632, 0.22294521665303071]
[0.7117642849034563, 0.18269016715718198, 0.45014059997829864, 0.37591414569832066, 0.24025174448410844, 0.29858693228235306, 0.23019800620416758, 0.44447201574787, 0.10749246293160268, 0.13870921713301645, 0.7156357593162649, 0.1551843231186634, 0.03740341134849236, 0.3985440860134391, 1.0161499732316146, 1.1156297864750935, 0.6804066357941937, 0.3959202823221901, 0.137583475539412, 0.2381829213237624]
[0.637864541406651, 0.1503421792257488, 0.3676657479332152, 0.320057566968519, 0.22564899934563057, 0.270500362707076, 0.17916060452846852, 0.33112405990653637, 0.10998421440894149, 0.14072106616155092, 0.6941158119081321, 0.15917588587457487, 0.03868647959581901, 0.40689888452777606, 1.0213288085687784, 1.1298756225554956, 0.7119244806213847, 0.4454799235433431, 0.14126170174114036, 0.25554938422424667]
[0.5559407848636867, 0.12631052480880015, 0.2904114839844158, 0.2582931568360568, 0.20220842330556968, 0.23220438765756146, 0.13519226638942444, 0.23330213347465842, 0.11226625726272813, 0.1400453134242365, 0.641570210214102, 0.1612004426742432, 0.03955097344654741, 0.40804286723593697, 1.0015444670165246, 1.1036887633479435, 0.7408738296753514, 0.5263705064271911, 0.14343585264469774, 0.2800916946303293]
[0.4860996729494896, 0.11179029473032814, 0.23137285784725128, 0.20552529077406026, 0.1747980995461694, 0.19357159968326387, 0.10612673454134991, 0.16719749268922784, 0.11358398642492225, 0.1352561059283022, 0.5613805721991025, 0.1598278497087544, 0.039596577742655764, 0.39819486451339087, 0.9479548560941534, 1.0301972322914708, 0.7517501261412793, 0.6414320506868165, 0.14297505072762562, 0.3074662732159045]
[0.43062336037116583, 0.10356398289144214, 0.1844615792588205, 0.16141813485406956, 0.14575039383099775, 0.1571093339964948, 0.08742077847876842, 0.12135214811929286, 0.11378070169199302, 0.1254693865013906, 0.4587858740179325, 0.15424893136862108, 0.03869286978329806, 0.37436167030100614, 0.8503781972412569, 0.9057553075538197, 0.729041383618697, 0.8111579252058726, 0.13952231909053012, 0.33452709185125584]
[0.3964192405936538, 0.10101479310886363, 0.1498840569184408, 0.12863185719239692, 0.12025445643257139, 0.1276930699613838, 0.07768609151302582, 0.09178438827929986, 0.1129872623811614, 0.11192393303962592, 0.3577855445251828, 0.14437936495554246, 0.03687729226745649, 0.3368923782528146, 0.7157846278632801, 0.7465856241369827, 0.6605821596120509, 1.0277574218259113, 0.133341718388981, 0.3523097808665832]
[0.38243072371733805, 0.10377513407300389, 0.1217115949038405, 0.10226707383008429, 0.09734334186196499, 0.1028555239811688, 0.07466666450788355, 0.07043649701997969, 0.11164051273424316, 0.09471976081247838, 0.2637113630658382, 0.12885618774692784, 0.033978371827748904, 0.28249894955757976, 0.544447847842159, 0.55639341927026, 0.53893465126863, 1.2848208126971767, 0.12395874464534769, 0.35570542529504545]
[0.3924286883443767, 0.11240697122613484, 0.10060160264185608, 0.08276595735749182, 0.07899615746867737, 0.08401465383098933, 0.0788422397090856, 0.05582188795693415, 0.11062078840452406, 0.07721940138253008, 0.19236961376758502, 0.10959561891383252, 0.03033101714965189, 0.22181251720504394, 0.38326331658163443, 0.3858474444456293, 0.40210575336666987, 1.478425292922194, 0.1124162514057935, 0.3416990564906213]
[0.42813571326082744, 0.128675511760532, 0.08395046251363436, 0.06718291407025544, 0.0632356350441606, 0.06877541295298652, 0.09176856327149403, 0.04453967283991596, 0.11135444226987187, 0.06048784203533888, 0.1407241774293773, 0.08818893948453095, 0.02603605400173823, 0.16352606263626918, 0.25543375688046666, 0.2550430395334262, 0.2841668822226807, 1.5316936451717071, 0.09862107795307937, 0.3119844172763782]
[0.4883402528108113, 0.15414578508205598, 0.07185284431714349, 0.05495095681805794, 0.04961797714957933, 0.05670337491344602, 0.11648832900242263, 0.0354977640435102, 0.11636390428998801, 0.0462790891202952, 0.10793925995233015, 0.06817276625365948, 0.021611980155946762, 0.11724461406725467, 0.1701205982863288, 0.17011786682827734, 0.20492293431105466, 1.4292032372179935, 0.0836268920392082, 0.2715005027005714]
⋮
[0.5101103239051296, 0.21403940327788046, 0.14189753714304945, 0.001693635023410558, 0.0, 0.001728937022754008, 0.28196217073426805, 0.0, 0.31327085058093246, 0.0, 0.09473818890539755, 0.0, 0.0, 0.0, 0.0004371189961953228, 0.16431048059903605, 0.4330100545383833, 0.16207440622211108, 0.0004404178929681443, 0.16833384448884556]
[0.5100179624486647, 0.21395516793778624, 0.142298340511395, 0.001399591351485802, 0.0, 0.0014285605990996871, 0.28201554962591546, 0.0, 0.3133389744030027, 0.0, 0.09462002403222947, 0.0, 0.0, 0.0, 0.00028712045202572195, 0.16437889628236, 0.4331528001886544, 0.1624194039181494, 0.0002894594976590179, 0.16827174480291335]
[0.5101391685673029, 0.21385402252043828, 0.14271524798887555, 0.001112636069065219, 0.0, 0.0011355054958041792, 0.28205082269333753, 0.0, 0.3133648590256169, 0.0, 0.09463567195834738, 0.0, 0.0, 0.0, 0.00017121921065168936, 0.16445747963247836, 0.4329436890213209, 0.16283785122750755, 0.00017269486463913314, 0.16827370452881724]
[0.5101982418484956, 0.21377046129833874, 0.1430518646451849, 0.0008782706978664726, 0.0, 0.0008962183476809098, 0.2820689756016431, 0.0, 0.3133576387193469, 0.0, 0.09471464863656844, 0.0, 0.0, 0.0, 9.947040358900451e-5, 0.1645130202313456, 0.43273362778296065, 0.1630429734495872, 0.00010035697721249949, 0.16835521456676708]
[0.5102084335756756, 0.2136967087621238, 0.143348139597704, 0.0006630399647620726, 0.0, 0.00067651590963556, 0.2820693689719039, 0.0, 0.31332800021060336, 0.0, 0.09477320176194699, 0.0, 0.0, 0.0, 5.144861940770425e-5, 0.16454897745245234, 0.4326529862755553, 0.16308117940893418, 5.191757472306722e-5, 0.1684581303519066]
[0.5101892000300914, 0.2136390751688376, 0.14358424821900018, 0.00048616264339761703, 0.0, 0.0004959989274277475, 0.28206673609400656, 0.0, 0.3133064900180514, 0.0, 0.09479046649120906, 0.0, 0.0, 0.0, 2.451430967046031e-5, 0.16456268461941692, 0.43267517418306095, 0.16307833641316483, 2.4740597400918844e-5, 0.16851125147059168]
[0.5101772405476588, 0.2135902678598142, 0.1437866089869059, 0.00033524651458835785, 0.0, 0.00034200273310671747, 0.28206907831200306, 0.0, 0.31329849477167543, 0.0, 0.09479928184557744, 0.0, 0.0, 0.0, 9.887997298776258e-6, 0.164565448997319, 0.43269238896886636, 0.16311370021270755, 9.979930181003285e-6, 0.1685218827160758]
[0.5101772405476588, 0.2135902678598142, 0.1437866089869059, 0.00033524651458835785, 0.0, 0.00034200273310671747, 0.28206907831200306, 0.0, 0.31329849477167543, 0.0, 0.09479928184557744, 0.0, 0.0, 0.0, 0.0, 0.164565448997319, 0.43269238896886636, 0.16311370021270755, 0.0, 0.1685218827160758]
[0.5101757757677363, 0.21355240023704158, 0.1439431400036674, 0.0002200025342075797, 0.0, 0.0002244227287650652, 0.2820723375375796, 0.0, 0.3132932906420358, 0.0, 0.09481482484946595, 0.0, 0.0, 0.0, 0.0, 0.16456927133861393, 0.4326792235117793, 0.16316037540338976, 0.0, 0.1685304597671773]
[0.5101703696591152, 0.213523666531707, 0.144060539077565, 0.00013284682506689382, 0.0, 0.00013550965212111025, 0.2820745504332237, 0.0, 0.3132897360866837, 0.0, 0.09482996255459025, 0.0, 0.0, 0.0, 0.0, 0.16456343881849694, 0.43267529562058404, 0.16317704041993888, 0.0, 0.16853030553320153]
[0.5101657341582638, 0.21350517347728484, 0.14413728087703948, 7.556735456477713e-5, 0.0, 7.707967503593065e-5, 0.28207544097209714, 0.0, 0.31328664374561793, 0.0, 0.09483755535192052, 0.0, 0.0, 0.0, 0.0, 0.16456072931990792, 0.43267635507709945, 0.16319154463635513, 0.0, 0.1685302162465592]
[0.5101635791665768, 0.2134928726933351, 0.14418784498203172, 3.80328307565163e-5, 0.0, 3.879316993321401e-5, 0.2820763883451882, 0.0, 0.31328518836452973, 0.0, 0.09484351749232077, 0.0, 0.0, 0.0, 0.0, 0.1645588376935647, 0.43267463374601806, 0.1632017694444765, 0.0, 0.16852982509285047]Step 3. Visualising the dynamics
Now we can move to plotting again. Note how we now ask for the time directly from the simulate object and all of the biomasses from that object as well.
Note too how we can supress the legend (which covers some of the time series).
plot(sim_niche, legend = false)
One game to play now is to alter the bodymass distribution. rand selects a randon uniform number between 0 and 1. Can you figure out how to make the distribution uniform between 0 and 10? See what that does.
A bit more about the process: dissecting the ModelParameters
Let’s dissect the ModelParameters object a bit, to understand just a bit more about what is going on.
params_nicheModelParameters{BioenergeticResponse}:
network: FoodWeb(S=20, L=80)
environment: Environment(K=[1, 1, ..., nothing, nothing], T=293.15K)
biorates: BioRates(d, r, x, y, e)
functional_response: BioenergeticResponse
producer_competition: ProducerCompetition((20, 20) matrix)
temperature_response: NoTemperatureResponseWalking through this
- The
networkcomponent defines the food web and reports the number of species and the links - the
environmentcomponent reports on values of the carrying capacity (K) and the baseline temperature (T). Note that K is specified only for the basal species via[1, 1, ... nothing, nothing]. All the producers have the same K at this point (1,1,1…). The presence of theTsuggests that we can ultimately work with climate change by running the model at different temperatures. There is a way to make some of the biorates and components of the functional response (see 3, 4) dependent not only on body mass, but also on temperature. - the
bioratescomponent contains detail on parameters central to making the model reflect a bit of reality:dis …;ris the intrinsic rate of increase (population growth rate) for the producers (plants);xandyare paramters associated with metabolic rates and consumption rates of the consumers (non-plant species). Finally,eis an efficiency/assimilation rate for the consumers eating either plants or other animals. - the
functional_responsecomponent defines the type of consumption function being used (e.g. Type I, II, or III sensu classic ecology and Holling). The functional response defines the interaction strength between species and how consumers change how much they eat dependent on the amount (density) of resource available.There are two options. The defaultBioenergetic Responsecorresponds to a phenomenological formulation where there are just two variables that describe how consumer consumption varies with resource density: a half-saturation parameter and an asymptote. The alternative calledClassic Responseis more trait based and includes the parameters attack rate and handling time. There are several other features of the functional response that can be manipulated, and these are introduced in later tutorials. - the
producer competitioncomponent allows us to alter how plants compete with each other. - the
temperature responsecomponewnt defines the absence or presence of temperature dependence, and when present, the shape of the relationship between biorates and functional response parameters and temperature.
Helper Functions: What can we do with the outputs?
As noted in the pre-amble, we are most often interested in additional information about the scenarios we build with the models. These include, for example, total biomass, biodiversity and stability. Let’s see how we can calcuate some of these. Keep in mind that to use these, we have to provide a model output AND a declaration of over how many time steps we want to calculate the metrics. For example, we probably want to calculate most of these metrics when the system is at equilibrium.
Before you start with the helper functions
In order to use these helper functions, we need to know the number of time steps that our model has produced.
It is important to know that ‘time’ is rather abstract in these models, partly because the integration of the equations to make the dynamics happens at a fine scale. In order to understand what the number of data points are in your simulation, you can look at the sim.t part of the outputs.
The other important thing to note is that the simulation have a built in stop point when the populations get to an equilibrium. This means that even if we were to ask for n= 100000 time steps, the simulations will stop when this ‘rule’ is met.
Finally, the algorithm is an adaptive one, so the number of output data points does not match ‘time’
We can find out what the number of steps is using the following code.
size(sim_niche.t) # this has 58 time points(58,)This little trick tells us the max number of time steps we can specify to calculate things, like the average biomass, stability, etc.
Using Some Helper Functions
First, we can get a meaure of total biomass in the community
# we use last = 30 because the max is 58
total_biomass(sim_niche, last = 30)2.487536153155013Second, we can an estimate of species persistence - how many have gone extinct! Remember that we started with 20, so a value of 0.6 means that there are 12 species left.
# here we specify last = 1 to get the estimate of % remaining from the start to end.
species_persistence(sim_niche, last = 1)0.6Third, we can look at stability, defined as the coefficient of variation. Here we look at using all of the data versus using data only out towards where the equilibrium is reached.
# here we can use ALL of the time steps
# size(sim_niche.t)[1] gets the actual number
population_stability(sim_niche, last = size(sim_niche.t)[1])-0.780564438364756Just for perspective, where you estimate stability is under your control. Here we estimate stability ‘out towards the equilibrium’ by specifying last = 30 rather than using the whole time series. This is the typical approach.
# here we use only the last 30.
population_stability(sim_niche, last = 30)-0.17360052844293952Finally, we can look at two measures of diversity. First, we can get species richness. This confirms that our persistence estimate (60%) is in line!
# con
foodweb_richness(sim_niche, last = 1)12.0And we can estimate evenness
foodweb_evenness(sim_niche; last = 10)0.8519930037484817What’s next
The next section will introduce how to create simulations where we work with multiple networks and collect data across these. We do this using loops, and collect information in a data frame. We then work on how to embed additional changes to parameters in the loops as well.